Melting of Polydisperse Hard Disks 
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The melting of a polydisperse hard disk system is investigated by Monte Carlo simulations in the semigrand 
canonical ensemble. This is done in the context of possible continuous melting by a dislocation unbinding 
mechanism, as an extension of the 2D hard disk melting problem. We find that while there is pronounced 
fractionation in polydispersity, the apparent density-polydispersity gap does not increase in width, contrary 
to 3D polydisperse hard spheres. The point where the Young's modulus is low enough for the dislocation 
unbinding to occur moves with the apparent melting point, but stays within the density gap, just like for the 
monodisperse hard disk system. Additionally, we find that throughout the accessible polydispersity range, the 
bound dislocation-pair concentration is high enough to affect the dislocation unbinding melting as predicted by 
Kosterlitz, Thouless, Halperin, Nelson and Young. 
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I. INTRODUCTION 

In the 1930's, Landau and Peierls showed that two dimen- 
sional solids are qualitatively different from their 3D coun- 
terpart, as they lack long-ranged positional order (see e.g. 
Ref. |1]). However, 2D crystals do have long-ranged bond- 
orientational order and, in this respect, they differ form the 
isotropic liquid phase where both translational and bond- 
orientational order are short ranged. 

In the 1970's, Kosterlitz and Thouless suggested that melt- 
ing of two-dimensional crystals may be quite different form 
3D melting. In particular, they proposed that melting in 
two dimensions may proceed via a continuous dislocation- 
unbinding transition. 

Kosterlitz and Thoulessj^, IH showed that the free energy 
associated with a single dislocation becomes negative when 
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where K is the Young's modulus of the crystal. At the point 
where K — IGnksT, dislocation pairs can unbind and, as solid 
with free dislocations 'flow' under shear, Kosterlitz and Thou- 
less interpreted this temperature as the melting point. In a 
more detailed analysis, Halperin and Nelson|2| and Young|4|, 
showed that dislocation unbinding is not enough to complete 
the melting process. At the point where the condition of Eq.^ 
is first satisfied, the system undergoes a (continuous) transi- 
tion from a 2D crystal to a hexatic phase. The hexatic phase 
is characterized by short-ranged (exponentially decaying) po- 
sitional order, but quasi-long-ranged (algebraically decaying) 
orientational order: the positional order is destroyed by the 
presence of the unbound dislocations. A second (continuous) 
phase transition is required to transform the hexatic phase into 
an isotropic liquid with short-ranged bond-orientational order. 

The Kosterlitz-Thouless-Halperin-Nelson- Young 

(KTHNY) theory makes precise predictions about the 
behavior of the coiTelation functions of both the translational 
and orientational order parameters. It should be stressed, 
however, that the KTHNY theory describes only a possible 
scenario: it also possible that one or both of the continuous 
transitions are first order, and even that there is a direct 
first-order transition from the crystal to the isotropic fluid. 



The KTHNY predictions sparked off an intensive search 
for real or model systems that would exhibit this two-stage 
melting process (for an early review see: ref. 15]. More recent 
examples can be found in refs. j^QI^ and Q). 

Surprisingly, however, there is still no satisfactory answer 
to the question whether the KTHNY scenario applies even to 
the simplest of all two-dimensional model systems, namely 
hard, elastic disks. In fact, this system was the very first to be 
studied in any computer simulation fioll . 

The reason why it is difficult to determine the nature of the 
melting transition, is that finite size effects tend to obscure the 
distinction between first order and continuous melting in 2D 
systems (see ref. f3\). 

In the case of hard disks, the early work by Alder and 
Wainwright suggested that the hard-disk melting transition 
was first order 1 11 1 (in the very early work of Metropolis et 
al. 1 10|, the computing power was insufficient to draw mean- 
ingful conclusions about the nature of the melting transition). 
The hard-disk melting problem was revisited many times af- 
ter the suggestion of the KTHNY scenario, but the evidence 
is still ambiguous. Evidence for continuous melting was re- 
ported in ref. 1121 . while evidence for a first-order phase tran- 
sition was presented in refs. fT3lfT4lfT5lfT6ll . In addition, sev- 
eral publications could not distinguish between the two sce- 
narios 1 17, 18, 19|. More recently, there has been some evi- 
dence for the KTHNY scenarioj^al^llQ but, flie matter still 
seems far from settled. 

One possible route to tackle this problem would be to con- 
sider hard disks as a special case of a more general class of 
systems, and study possible trends in the phase behavior of 
this generalized model. In the present case, we consider the 
2D hard-disk system as a special case of polydisperse disks. In 
3D hard spheres, the melting transition is of first order. As the 
polydispersity is increased the difference in volume fraction 
of the coexisting solid and liquid phases widens with increas- 
ing polydispersity 1 23 , 24,1 . One might expect similar behavior 
in two dimensions if the solid-liquid transition would be of 
first order. 

While polydispersity in two dimensional systems has been 
studied before, it was in the context of melting by increasing 
size dispersity for Lennard Jones systems Ii25il or in the con- 
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text of a possible glass transitionf26, 27]. In the present paper, 
we examine the phase behavior of polydisperse hard disk sys- 
tems. 



II. THE SYSTEM 

A. The Semigrand Canonical Ensemble 

The model for polydispersity is based on the semigrand 
canonical ensemble|281; this ensemble has previously been 
used to study the phase diagram of polydisperse 3D hard 
spheres! 23, 24]. The semigrand canonical ensemble can be 
seen as a hybrid version of the canonical ensemble and the 
grand canonical ensemble. It is characterized by a thermo- 
dynamic potential X that satisfies the following fundamental 
thermodynamic relation L29il 



dX = -SdT-PdV- J N{a)?>p{a)da 



(2) 



Here S is the entropy of the system, T the temperature, P the 
pressure, V the volume. N{a)da denotes the number of par- 
ticles with diameter between a and a + da and /J (a) is the 
chemical of particles with diameter a. We now add and sub- 
tract a term containing the chemical potential of a reference 
species Oq from the complete differential: 

dX = -SdT - PdV-Ndij{ao) - [ N{a)5Aij{a)da (3) 



where we replace /j(o) — ^((oo) with A/j(a). We now perform 
a Legendre transformation to a new ensemble that has as 
a thermodynamic control parameter instead of /j(oo) (and P 
instead of V). 



dY = -SdT + VdP + ij{aQ)dN- J N{a)dAij{a)da 

which, in explicit form becomes (with the Euler equation) 

Y{N,Aij{a),P,T) = U -TS + PV+Nij{ao) 
-J N{a)Aij{a)da 
= N/jioo) 

The partition sum for this ensemble is 



(4) 



(5) 



Y{N,Afi{a),P,T) 



exp (-P {PV + U{V, s'^) - A/j{a)N{a)]}) 



da^ dV ds 



(6) 



with Y = —ksTlnT. For simulation purposes, the semigrand 
canonical ensemble can be interpreted as one where there is 
constant number of particles that can change identity. This 
identity switching can be an extra move in a Monte Carlo 
simulation; in a polydisperse mixture, this would amount to 
a particle size change with an acceptance criterion based on 
the functional form for the chemical potential A/j(a). 



As in Ref. 112 3r . we use the following functional form for 
the chemical potential 



Mo) 



-k„T- 



2v2 



(7) 



which, at zero density, will give a Gaussian particle size distri- 
bution according to the partition sum of Eq.|S] In practice, the 
size distribution is Gaussian-like at the densities of the crys- 
talline phase. 



B. Order Parameters 

For the study of the properties of the phases and the phase 
transitions, some order parameters have been used that are 
standard in the studying of 2D melting|5| and for which the 
KTHNY melting scenario makes explicit predictions! 2|. We 
define the n-fold bond-orientational order at x,, the position of 
particle /, as 



(8) 



where A^, is the number of neighbors and 0; (x') is the angle 
between an arbitrary (fixed) axis and the line connecting par- 
ticle / with its y-th neighbor; two particles are neighbors if they 
share a Voronoi cell edge. For systems that tend to crystallize 
into triangular lattices, the leading bond-order parameter is the 
one for which n = 6. The global value of the order parameter 
is simply the mean of the local values. 

The positional order is measured using the static structure 
factor S{q) at one specific scattering vector q equal to a re- 
ciprocal lattice vector of a perfect crystal with orientation and 
lattice spacing taken from the system. To check for hexagonal 
crystalline positional order, the lattice vector oq is set to its 
ideal value for a given packing fraction: 
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The crystal orientation is taken from the mean angle obtained 
from the global hexagonal bond-orientational order parame- 
ter: 



(10) 



which specifies the orientation of one of the six equivalent 
crystal axes within an angular range < a < 7i/3. Once the 
average orientation of the nearest-neighbor 'bonds' ao has 
been specified, it is straightforward to deduce the orienta- 
tion of the corresponding reciprocal lattice vector G through 
G • ao = 2k. The positional order parameter ^ of the crystal is 
then given by 



C(x,) 



JGxj 



(11) 
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Radial correlation functions of the order parameters are de- 
fined as 



C(r) = (r(0)C(r))/g(r) 



(12) 
(13) 



In two-dimensional systems, ^ is expected to decay to zero, 
either exponentially (short ranged order) or algebraically 
('quasi-long ranged'). The KTHNY melting scenario makes 
predictions for the decay of the orientational order in the hex- 
atic and phase: gei^) ^ r^^* with Tjg j at the melting of 
the hexatic phase into the liquid phase|0l3. 



C. Elasticity 

Eq.^ provides a very useful test to decide whether or not 
a 2D melting can be of the KTHNY type. If we find that 
Young's modulus drops below the 'magical' value of 1671^^7 
in an otherwise stable solid, then it is very likely that this solid 
melts by dislocation unbinding. Conversely, if we find that 
this magical value is only crossed at densities where we know 
that the isotropic liquid phase is thermodynamically stable, 
then it is reasonable to assume that melting is a first-order 
transition. Often, however, the simulations do not provide a 
clear answer, as the point where K = 1671^^7 is located in the 
intermediate density regime that may either be a two-phase re- 
gion separating two stable phase, or the domain of the elusive 
hexatic phase. 

The Young's modulus is defined through the shear (A-l) and 
bulk {pi} Lame elastic constants in 2D, 



K = 



2^L + Xl 



< 16% 



(14) 



where ciq is the equilibrium lattice spacing. The Lame elastic 
constants are related to the second-order elastic constants (i.e. 
the elastic constants defined by the second derivative of the 
free energy to the Lagrangian strain), through 



Cii = Xl + ^/jl 
Cn = Xl 
Caa = ial-P 



(15) 



The theory of the KTHNY scenario is based on a group 
renormalization of both the Young's modulus and the disloca- 
tion fugacity 



y = exp 



(16) 



where Ec is the core energy of a dislocation. As the length 
scale of the renomalization is increased towards infinity, the 
dislocation fugacity and the elastic constants that constitute 
K tend to their 'renormalized' values, where the (long-range) 
effect of (presumably low-concentration) dislocations is prop- 
erly taken into account. 

The elastic constants measured in this work, although cal- 
culated at finite sizes, should be close to their renormalized 



values because of the high concentration of dislocations in the 
hard disk system, as will be shown in Sectiorlllll 

When specifying elastic constants of a polydisperse system, 
we should distinguish between 'quenched' and 'annealed' 
elastic constants. Quenched elastic constants measure the sec- 
ond strain derivative of the free energy of a polydisperse crys- 
tal with a 'frozen' size distribution: i.e. the particle-size dis- 
tribution is assumed not to respond to the deformation. In 
contrast, the 'annealed' elastic constants measure the second 
strain derivative of the semi-grand potential. In this case, the 
particle size distribution is assumed to respond to the defor- 
mation. The quenched constants are the ones that are pre- 
sumably measured in mechanical experiments that probe the 
elastic deformation of a polydisperse solid. But the annealed 
constants describe the equilibrium state of a deformed poly- 
disperse solid. In order to determine the critical value for 
Young's modulus in a polydisperse 2D solid, we computed the 
annealed elastic constants, as these determine the equilibrium 
behavior of the system. 

To calculate the elastic constants at different polydispersi- 
ties, a hybrid Monte Carlo - Molecular Dynamics method was 
used, where Monte Carlo runs sampling particle diameters 
and positions in the strained system were followed by Molecu- 
lar Dynamics runs where the stress tensor was measured using 
the instantaneous realization of the particle size distribution. 
The algorithm thus effectively calculates the elastic constants 
of many realizations of a polydisperse hard-sphere crystal. 

The quenched elastic constants will be larger than the an- 
nealed quantities, because of the concavity of the free energy. 
Hence the corresponding 'quenched' Young's modulus will 
only reach the instability limit K = 16k at lower densities. 

The strains used in the elastic constant determination were 



a^ = 



l+a 


1 b 
1 





l+a 



(17) 



with a as the (small) strain parameter The resulting stress 
derivatives give us tSOfl 



dTn 
da 

dTn 
db 



dTu ar, 
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da 



= Cu + Ci2 = 2Xl + 2^L 
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-Caa+P-- 



(18) 
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allowing us to calculate the Young's modulus of Ea.[T4l 



III. SIMULATIONS 

The simulations were performed on systems containing 
59 X 68 = 4012 particles, with a chemical potential disitri- 
bution width ('polydispersity') parameter v (see Eq.Q vary- 
ing between v — 0.00025 and v — 0.008. For higher values 
of V (higher polydispersities), the equilibration was exceed- 
ingly slow, even with the combined volume-particle radius 
sampHng: equilibration took 1.3 x 10^ to 2.5 x 10^ MC steps 
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FIG. 1 : Calculated equations of state (with normalized pressure) for 
varying polydispersity control parameter v. System size is 4012 par- 
ticles. 

per particle and data were sampled during 1x10^ steps per 
particle; around 20% of steps were devoted to particle radius 
steps. For = 4012, the simulation time was 6-8 hours on an 
Athlon 1600+ CPU. 

Monte Carlo simulations of the semigrand canonical en- 
semble of Eq.|6] involve displacement moves, volume moves 
and identity change moves, which in the current context mean 
particle radius moves. To speed up the simulation, the cou- 
pling between mean particle radius and system size is ex- 
ploited to remove the global volume moves and simultaneous 
systemH-particle size scaling was introduced. This improves 
the statistics of sampling considerably! 23 1. It turns out that, 
in practice, a real simultaneous system-nparticle size scaling is 
at least as fast as the analytical integration of the scaling part 
of the partition sum also described in Ref. |23]. 

The resulting equations of state are shown in Fig. [2 
From this figure, it is clear that the phase transition, which 
is around packing fraction r| w 0.70 for the monodisperse 
easel ITlllSlflQll . shifts to higher packing fractions and higher 
pressures with increasing polydispersity. At v = 0.008, the 
system cannot be made to freeze at all. A similar phenomenon 
was observed in ref. |23] — it is due to the choice of the func- 
tional form for the chemical potential (Eqn. 0. In addition, 
for V = 0.007, it is extremely difficult to equilibrate the sys- 
tem properly in the vicinity of the phase transition. 

As is the case for 3D spheres, the liquid branches of the 
equations of state very nearly superimpose i23ll . However, 
near the phase transition, the pressure appears to increase 
slightly with polydispersity. Upon further compression, the 
system undergoes a phase transition with an apparent density 
gap. However, because of the relatively small system size, 
the presence of such an apparent density gap is also compat- 
ible with continuous melting. We find that the density gap 
decreases with increasing polydispersity. 

In Fig. 12 we plot the variation in polydispersity upon 
freezing. As a measure for the polydispersity, we use s = 
(o2)/(o>2- 1. In the same figure, the apparent density gap is 




FIG. 2: Packing fraction r] as a function of polydispersity s for vary- 
ing V, with W = 4012. The boundaries of the apparent density gap are 
shown as the solid lines; the density gap of the monodisperse (s = 0) 
case is taken from data by Jaster| 18], with N = 4096. The dashed line 
shows the location of packing fractions (with matching polydispersi- 
ties) extrapolated to be A" = 16lt; the gray area denotes an estimate 
of the size of the error. The black circle shows the location of the 
A" = 1671 point obtained using a larger system size (A' = 4012) for 
the elasticity simulations. 




FIG. 3: Mean particle diameter (a) as a function of packing fraction 
r| for the different values of the polydispersity control parameter v. 



also shown (the borders of the density gap shown here are sim- 
ply the solid with the lowest packing fraction and the liquid 
with the highest pack ing fraction found in the simulations). 
As in the 3D case I23ll24ll . the freezing point moves to higher 
volume fractions as the polydispersity is increased and size 
fractionation is also increased. But, whereas the density gap 
widens upon freezing in 3D, it appears, if anything, to de- 
crease in 2D (except for the possibly seemingly poorly equi- 
librated case of V = 0.007). Additionally, the maximum poly- 
dispersity at which the solid seems to remain stable (around 
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FIG. 4: Example of the decay of the orientational (left) and positional 
(right) order parameters over the phase transition, for v = 0.006. The 
decay of the positional order parameter, ^6('') goss as gs(r) o= r^4 at 
the intermediate pressure -P(CT^) = 12.7. The positional order param- 
eters decay exponentially for both lower pressures. 



8%) is considerably higher than in 3D (5.7%). 

In the semigrand-canonical ensemble, the system equili- 
brates to a mean particle size, dependent on pressure and v. 
Fig- E shows this mean particle size as a function of packing 
fraction: the phase transitions show up as jumps in the pack- 
ing fraction, but not in the mean particle size. This means that 
there is no particle size fractionation (the mean of the particle 
diameter distribution does not change) while there is, as can 
be seen in Fig|2l polydispersity-fractionation (the width of the 
particle diameter distribution changes). 

An example of the behavior of the order parameters near 
the phase transition is shown in Fig.|3 The positional and ori- 
entational order both increase sharply (but not quite simulta- 
neously) in the region of the phase transition, similar to what 
one finds for monodisperse hard disks. Although the decay 

of the orientational correlation function g(,{r) goes with r^? 
as predicted by the KTHNY scenario for the hexatic at the 
hexatic-crystal transition, this seems to be a coincidence, be- 
cause a further simulation of the same configuration yields a 
different (but still algebraic) decay rate. This indicates that the 
typical decay 'time' of fluctuations in the system exceeds the 
typical length of the (rather long) simulations. 

We stress that all observations reported thus far are con- 
sistent with either a KTHNY scenario or a weak first-order 
transition. The KTHNY scenario however, is based on the as- 
sumption that the concentration of bound dislocations in the 
crystal phase is low; dislocation interaction is not taken into 
account and the Kosterlitz-Thouless normalization is based on 
an expansion in the dislocation unbinding length which may 
be unrealistically long for high dislocation concentrations. We 
measured the concentration of bound dislocation pairs (see 
Fig. In the figure, we show the concentration of seven- 
coordinated particles. Because in the crystal the number of 
eight or more coordinated particles is negligible, and as the 
number of point defects turns out to be an order of magni- 
tude lower than the number of bound dislocations^!^, the 
number of seven-coordinated particles is a good measure of 
the dislocation count. At the melting point (the boundary of 
the apparent density gap) the dislocation concentration varies 
from 1 % to more than 3%. As the concentration of dislocation 
pairs depends sensitively on the dislocation core energy, this 
suggests that the core energy is rather low. Note that a core 
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FIG. 5: Concentration of seven-coordinated particles xj (a measure 
for the number of dislocations) as a function of packing fraction rel- 
ative to the melting packing fraction (determined by the boundaries 
of the apparent density gap). 



energy less than 2 to 4-kBT is not compatible with KTHNY 
meltingiHlll. 

Neglecting the dislocation-dislocation interaction (a very 
crude approxximation), but including the internal elastic en- 
ergy and entropy of dislocation pairs|20, 33], this defect den- 
sity would be consistent with a core energy of approximately 
SksT over the full polydispersity range. This, however, is a 
crude approximation and the real core energy may be several 
kgT off. The value is compatible with both the value of GksT 
given in Ref |20] and the value of llksT at a much higher 
packing fraction of 0.82 calculated in Ref. i3 . 



IV. ELASTIC CONSTANTS 

The elastic constants were measured using the method de- 
scribed in section III CI using simulations of 412 particles, 
equilibrating for 4 x 10'' MC steps and 4 x 10'* MD colli- 
sions per particle and measuring up 5 x 10^ MC steps and 
1.5 X 10^ MD collisions per particle. Earlier work 1 19] had 
shown that the elastic constants are not significantly affected 
by the presence of point defects such as vacancies. Away from 
the KTHNY transition, the elastic constants are not very sen- 
sitive to finite-size effects |14]. Of course, this is not true 
close to a KTHNY transition, but this will turn out to be less 
relevant for the present system because we always observe 
melting before we get into the 'danger zone' . 

The results of the simulations are shown in Fig. |5J here 
an exponential of the Young's modulus, K = (1 — 16k/ Ky*"^, 
where Vkt is the Kosterlitz-Thouless renormalization expo- 
nent, which has a value of Vkt ~ 0.3696312]. This form is 
chosen because close to the KTHNY transition, the Young's 
modulus — expressed as function of the difference between 
the packing fraction r| and the packing fraction at the KTHNY, 
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FIG. 6: Fractional exponent of the Young's modulus K = (1 — 
167i/A:)^kt with vk^t = 0.3696 () for the different polydispersities, 
as a function of packing fraction. The lines are the fits to determine 
the locations of A" = 16jr. The points at v = 0.004 with the trian- 
gles together with the fit shown with the dashed line are results of 
calculations done with 4012 particles. 

r|KT — behaves asjH 



1671 1 -c(r|-riKTrKT 

The measured K values are then fitted to a second-order poly- 
nomial. The results for the monodisperse case (v = 0) are 
similar to those by Wojciechowski and Branka fl4l Bsll and 
by Bates and Frenkel L19tl . 

To check for finite-size effects in the elastic constants cal- 



culations, the calculations for v = 0.004 were also done with 
4012 particles. The results are shown in Fig.|6land in Fig.|2l 
the = 167t extrapolation is very similar to the smaller 412 
particle system and agrees within the error margin. 

From the locations of the K — 167t line in Fig. Hit is clear 
that the situation with respect to the type of phase transi- 
tion is, even for higher polydispersities, similar to that of the 
monodisperse system; the points seem to follow — within the 
statistical uncertainties — not only the loci of the phase tran- 
sitions, but also the position within the density gap. 



V. CONCLUSION 

In this article, we explored the effect of polydispersity on 
the nature of the 2D melting transition is a system of 2D hard 
disks. 

We find that the solid-liquid phase transition shifts to higher 
packing fractions as the polydispersity increases, and that 
polydispersity fractionation takes place in the region of the 
phase transition. The maximum polydispersity at which the 
solid can be stable is larger than in 3D hard spheres. 

The density-polydispersity gap, be it real or apparent, does 
not seem to increase in size with increasing polydispersity. 
The fact that the points for which K —Xdn appear to be lo- 
cated in the two-phase region, supports the assumption that 
the melting transition is first order Even if this should not be 
the case, the high dislocation concentration will presumably 
have an effect on the KTHNY predictions. 

The work of the FOM institute is part of the research pro- 
gram of the Foundation for Fundamental Research on Matter 
(FOM) and was made possible through financial support by 
the Dutch Foundation for Scientific Research (NWO). 
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